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SOLUTION  OF  THE  PERTURBED  HARMONIC 
OSCILLATOR  BY  COMPUTER 

I.  INTRODUCTION 

Many  problems  in  applied  mathematics  reduce  to  a non- 
linear oscillator.  In  fact,  in  celestial  mechanics  many 
times  transformations  are  sought  which  will  reduce  a prob- 
lem to  a nonlinear  oscillator.  Some  examples  are: 

1.  The  central  force  problem  after  the  transformation 

2 

r * (1/u),  dt  = r dr  where  r is  the  radius  and  t is 
t ime ; 

2.  Motion  in  the  neighborhood  of  a stable  equilibrium 
point  and 

3.  The  Vinti  problem  [1],  motion  about  an  oblate 
spheroid,  which  can  be  reduced  to  two  uncoupled  perturbed 
harmonic  oscillators  and  a quadrature.  The  purpose  of  this 
report  is  the  development  of  an  analytic  series  solution  of 
the  perturbed  harmonic  oscillator 

*X  + X + €f(X;£)  - 0 (1) 

where  f(X)  is  a polynomial  in  X and  £ is  a small  para- 
meter . 


Not«;  Manuicript  tubmitted  April  28,  1977 
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Low  order  solutions  can  be  obtained  "by  hand"  but  the 
amount  of  algebraic  manipulation  required  to  obtain  a high 
order  solution  makes  carrying  the  solutions  to  higher  order 
prohibitive.  This  report  describes  a computer  program 
that  develops  an  analytic  series  solution  to  the  perturbed 
harmonic  oscillator,  Eq.  (1),  to  the  desired  order  in  e. 

Lindstedt's  method  [2]  is  used  to  develop  the  solution 
and  algebraic  manipulation  routines  13]  are  employed  to 
carry  out  the  algebraic  manipulation. 

II.  TECHNICAL  DISCUSSION 

A.  Motion  Near  an  Equilibrium  Point 

As  an  example  of  a problem  which  reduces  to  a 
perturbed  harmonic  oscillator,  consider  the  problem 

y = F(y)  , F(y^)  « 0 , |^  <0  (2) 

y-7o 

Thus  y * y is  a stable  equilibrium  point.  To  analyze  the 
o 

motion  in  the  neighborhood  of  the  equilibrium  point,  expand 
F(y)  in  a Taylor  series  about  y^,  i.e, 

F(y)  - F(y^)  + F'(y^)(y-y^)  + jF ' ' (y^) (y-y^) 2+  ...  (3) 

and  let 

X - (y-y^)/L  (4) 

where  L is  a representative  scale  length  so  that  x * 0(1) 
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at  the  maximum  value  of  X.  With 
e = 1/L 

(5) 

(7^  * - 

Equation  (2)  reduces  to 


X*  + a^X  * - ef(X;e)  (6) 

where  f(X;e)  is  polynomial  in  X.  A further  time  transform- 
ation of  s = 7t  transforms  the  problem  to  one  with  a 
frequency  of  unity. 

An  exact  solution  of  Eq,  (1)  can  be  obtained  in  terms 
of  elliptic  functions  but  sometimes  elliptic  function  solu- 
tions are  not  particularly  desirable  and  the  reduction  of 
the  problem  to  obtain  the  elliptic  function  solution  can  be 
quite  difficult.  There  are  several  methods  for  obtaining  a 
series  solution  in  power  of  e but  the  best  methods  for 
this  problem  are  Lindstedt's  method  [2]  and  Lie  Series  . 

Lie  Series  involves  the  use  of  Hamiltonian  mechanics 
and  is  easily  automated  on  the  computer.  For  studying 
properties  of  the  solution,  it  is  probably  the  best  method, 
but  to  obtain  solutions  in  terms  of  specific  initial  condi- 
tions, Lindstedt's  method  is  the  easiest  to  implement. 

B.  Lindstedt's  Method 

Consider  the  problem: 
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(7) 


X + X = - ef(X;-:),  X(0)  = A,  X(0)  = B 
To  develop  a series  solution  let 


x=i: 

n=o 


n 


X = 
n 


e^X^  + 


(8) 


Since  the  frequency  is  amplitude  dependent  in  nonlinear 
oscillator  problems,  introduce  a new  time  scale  t via 


= JUt 


(1  + 


2 

+ e -*■'2 


)t 


(9) 


where  the  x'^  are  undetermined  at  this  point.  Substituting 
Eqs.  (8)  and  (9)  into  Eq.  (7)  and  equating  like  powers  of 
c gives 


X" 

4. 

X = 

0 

0 

0 

Xf 

+ 

II 

X 

-2ju,X"  -f(X  ;0) 
1 0 0 

/ / 

X 

n 

+ 

X - 
n 

®n^'“l’'“2’ • • • ’ Vl’ 

X 

n- 


1 


) 


(10) 


The  initial  conditions  are  satisfied  by  the  zeroth  order, 

i.e. 


X^(0)  = A,  X^(0)  - 0,  n » 1,  2,  ... 
X^(0)  - B 

X{(0)  + T^X'^(O)  - 0 

X'(0)  + -x),X'  .(0)  + ...  + o)  ,X'(0)  - 0 
n 1 n- Jl  n- i o 

4 


(11) 


(12) 


J 


Eqs.  (12)  are  a result  of  setting  the  zeroth  order  of 
" = equal  to  B and  the  higher  orders  in  € equal 

to  zero. 

The  solution  of  X is 

o 

X = Acost  + Bsin*  (13) 

Substitution  of  Eq,  (13)  into  the  X^^  equation  results  in 
resonance  of  the  X^  equation  because  the  homogeneous  portion 
of  Eqs.  (10)  are  identical.  Since  a uniformly  valid  solu- 
tion is  sought,  no  resonant  or  secular  terms  can  be  present 
inX,n=l,  2,  ....  Therefore,  d,  is  chosen  to  make  the 
coefficient  of  the  resonant  term  in  X^  zero.  The  rest  of 

the  D terms  are  evaluated  in  like  manner, 
n 

C . Example 

As  an  example,  consider  a restricted  form  of  the 
Duffing  equation,  i.e.,  the  hard  spring  problem 

X + X + eX^  = 0,  X(0)  = A,  X(0)  = 0 (14) 


Application  of  Eqs.  (8),  (9)  and  (12)  yields 

ef 

X"  + X = 0, 
o o ’ 

X (0)  = A,  X (0)  - 0 

0 o 

3 

X."  + X,  - -2m;,  X"  - X 

1 1 loo 


(15) 


(16) 
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Xj^(O)  = X{(0)  = 0 

x''  + X„  = -(2X'  + xh  x"  - 2u),x,"-  3X^  X. 

2 2 2 l o 11  ol 

X2(0)  = X^  (0)  = 0 


(17) 


The  solution  to  the  zeroth  order  equation  is 

X = Acost  (18) 

o ^ 

Substitution  into  Eq.  (16)  gives 

3 

X/ ^ + X^  = A(2x-  - y A^)cos-^  - ^ cos3^ 

11  14  4 

The  cos~  term  results  in  resonant  or  secular  terms  in  X^ 

which  cannot  be  present  if  a uniformly  valid  solution  is  to 

be  obtained.  To  suppress  this  secular  term  set 


^ A^ 
8 ^ 


(19) 


The  corresponding  solution  of  X^  which  satisfies  the 
initial  condition  is 

a3 

Xj^  = -jj  (cos3t  - cost) 

The  solution  of  Eq.  (14)  is 

13  0 

X - Acos"  + (cos3’’  - cost)  + 0(€^) 


1+1  cA^  + 0(e2)~  t 

O J 


(20) 


(21) 

(22) 
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III.  COMPUTER  IMPLEMENTATION 


The  computer  program  LNDSTD  has  been  written  to  carry 
out  the  Lindstedt  procedure.  The  program  currently  runs  on 
the  NRL  DEC-System  10  computer  and  is  available  on  request 
from  Code  7904.  Inputs  to  the  program  and  other  considera- 
tions are  given  below.  Examples  are  given  at  the  end. 

A.  Algebraic  Expressions  Written  as  Poisson  Series 
Before  specifying  the  inputs  to  LNDSTD,  it  is 
necessary  to  describe  the  way  in  which  algebraic  expressions 
are  coded  in  the  algebraic  routines  used.  LNDSTD  uses  the 
routines  of  R,  Dasenbrock  which  are  described  in  detail  in 
ref.  (3).  In  order  to  use  these  routines,  an  algebraic 
expression  must  be  written  as  a Poisson  series  in  which 
each  term  is  of  the  form; 


■ 12 

' 6 

N 

D 

TT 

\cos  J 

y a. 9. 

jL^  1 1 

i-l 

_i=l 

where : 

N and  D are  decimal  integers  up  to  11  digits; 
v^  are  polynomial  variables; 
e^  are  the  corresponding  integer  powers; 
are  trigonometric  variables; 
are  integer  coefficients. 
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Each  term  is  coded  using  six  computer  words.  The  first  two 
words  contain  N and  D respectively  entered  as  decimal 
integers.  The  exponents  e^  - e^^  corresponding  to  variables 
v^  - v^2  3re  packed  into  words  three  and  four  as  two  octal 
integers.  (Note  that  an  octal  integer  is  indicated  both 
here  and  on  the  DEC  System  10  by  using  double  quotation 
marks  preceding  the  number) . The  fifth  word  is  +1  if  the 
cosine  function  is  intended  and  -1  if  the  sine  function  is 
intended.  The  coefficients  are  packed  into  the 

sixth  word  as  an  octal  integer.  The  way  in  which  a term  is 
coded  into  these  six  words  is  best  explained  by  the  follow- 
ing example: 

3 3 -2  9 11 

Code  the  expression  -r  v,  v,  v.  v„  cos  (9,  + 59,) 

41458  1 6 

into  six  computer  words.  The  expression  makes  up  a Poisson 

series  of  one  term  and  may  be  written  as 

3300  -2  900  11  00  0 0 

4 '^l  ^^2  "^3  ^4  ^5  ^6  '^7  ^8  "^9  "^10  ^11  "^12 

. X cos  '’0,  + 0e„  + 09-  + 09,  + 09.  + 59. 

-i  2 J 4 5 D. 

Word  1 is  set  equal  to  3. 

Word  2 is  set  equal  to  4. 

Each  two  octal  digits  of  words  3,  4 and  6 are  set  equal  to 

the  octal  value  of  the  correspondong  exponents  or  coeffi- 
cients biased  by  the  amount  ”40  for  variables  ~ , 
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^1  ” ^2  ” ^6  amount  ”20  for  variables  v^ , 


Word  5 is  set  equal  to  +1. 
The  result  is: 


Word  Word  Word 

12  3 

Word 

4 

Word 

5 

Word 

6 

N D 

®1®2®3®4®5®6 

^7®8®9®nfll®L2 

sin 

cos 

^l'^2^3^4‘^5^6 

3 4 

”234040365140 

”205340404040 

1 

”214040404045 

In  what 

follows,  any  reference  to  a variable  v 

. or  S . is  an 

1.  L 

implicit 

reference  to  the 

position  of  the  corresponding  e^ 

or  a.  within  word  3,  4 or 
1 ’ 

6. 

B. 

Inputs  to  LNDSTD 

1.  The  first  step  in  running 

LNDSTD 

is  to  specify 

whether  input  will  be  from  the  terminal  or  from  disk  and 
whether  output  will  be  to  the  terminal  or  to  disk.  In 
both  cases,  the  response  'T'  indicates  the  terminal,  and 
any  other  response  indicates  disk.  If  disk  is  chosen  in 
either  case,  a file  name  will  be  requested.  Note  that  the 
extension  .DAT  is  assumed  for  the  file  unless  some  other 
extension  is  specified  or  the  file  name  is  ended  with  a 
period. 

The  input  file  should  be  prepared  following  the 
formats  of  steps  2-7  below.  All  line  sequence  numbers  must 
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be  removed  from  the  file. 


2.  The  next  input  record  specifies  the  desired 

order  in  of  the  solution.  The  order  may  not  exceed  4. 

3.  The  next  input  record  specifies  the  number  of 

terms  in  f(x;€),  (NTERM) . 

4.  The  next  10  input  records  specify  names  for  €, 
X,  initial  value  of  x,  initial  value  of  x,  and  5 other 
variables  available  to  the  user  for  defining  f(x;€). 

Carriage  returns  may  be  entered  for  variables  which  will 
not  be  used.  Names  may  be  up  to  four  characters  long. 

5.  The  next  NTERM  records  define  f(x;-:).  Note 
that  the  f(x;€)  used  must  correspond  to  the  differential 
equation  written  in  the  form 

x + x+-cf(x;5)=0  (24) 

Each  term  is  entered  as  one  record.  The  function  f(x;-:) 
must  be  a polynomial  in  x and  may  not  contain  any  trigo- 
nometric variables,  i.e.,  each  term  is  of  the  form 


N 

D 


(25) 


Thus,  each  term  is  composed  of  a numerator,  a denominator, 
and  two  octal  words  corresponding  to  variables  1 to  12.  If 
c appears  in  a term,  variable  Vj^  must  be  used  for  it. 
Variables  v„  - v,  must  not  be  used;  i.e.,  they  should  be 


10 


,,  be  used  £ot 

^ a<?  ”^0's.  Varxa  12  format 

entered  - v,,. 

•oKleS  should  D 

other  varrable  _ 

h term  i-S  (2U0.  ’ of  the 

jor  each  tern  ^be  fc^ 

f.  The  next  r 

in  general,  the  form 

3,no  order  solutrcn.  * 

^ i <;T.n  'S 


^ ^ V S lU  ^ 

7rl  the  initial  values  o£  h -d  x 
where  x^  and  x^  solution 

respectively,  ^nter  ^e  sine  solution  is  desire 

,.  _ OV  'S'  ° . is  desired 

desired  (^i  ' ’ ine  solutron 

1 • <:  c»  srne  anu  v- 

(x.^  = 0)-,  or  'SC  t 

(Xb  . 0 and  X,  - 0)  . ^ „^rions  avail- 

The  last  inpu  ^re  input 

able  to  the  user.  Six  option^  ^ ,3 

tn  the  £or.at  («)•  ^nter  ^ 

,esired  and  aero  otherwise. 

bent  to  choosins  no  options  ^ „^,bon 

1)  Skip  the  c ec  substitute  the 

the  progran  writ 

is  not  chosen,  ^.^ge.ential  equation 

final  solution  into 

. If  the  solutxo 

as  a check.  executxon 

re  .ore  store. 


Unal  check  ®ay  storage 

oosslbly  heduira  ">nT^ 
tii„e  and  may  P°ss 
1 ^ 


u 


2) 

Do 

not 

type 

final  X. 

3) 

Do 

not 

type 

final  X. 

4) 

Do 

not 

type 

final  '. 

5) 

Do 

not 

type 

intermediate  results. 

6) 

Decide 

at  each  step  whether  to  type 

mediate  results.  If  this  option  is  chosen, 
LNDSTD  will  type  out  at  each  output  step  the 
number  of  terras  in  the  expression  corresponding 
to  that  step  and  will  then  type  ”PRINT(Y/N) ?" 
The  user  should  type  "Y"  to  have  the  expres- 
sion typed  out.  Any  other  response  will 
suppress  typing,  and  processing  will  continue. 
If  the  first  number  typed  is  other  than  a 0 or  a 1,  this 
list  of  option  descriptions  will  be  typed. 

C . Output  from  LNDSTD 

LNDSTD  first  types  out  the  information  input  by 
the  user  so  that  it  may  be  verified.  (The  user  is  asked 
whether  or  not  to  proceed.  A negative  response  causes 
LNDSTD  to  start  over  from  the  beginning.)  If  none  of 
options  2-6  were  taken,  the  program  proceeds  to  type  all 
intermediate  and  final  results.  Intermediate  results  con- 
sist of  the  solutions  at  each  order  i of  t'  and  x.. 

i L 

Final  results  consist  of  the  solutions  of  x,  x,  and  + to 
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the  order  in  specified  in  step  2 of  the  input.  If  any 

of  the  options  2-6  are  chose,  output  of  results  to  the 
terminal  is  changed  accordingly.  If  option  1 was  not 
chosen,  LNDSTD  will  carry  out  a check  on  the  solution  by 
evaluating  the  expression  x + x -f-  -:f(x).  If  the  check  is 
successful,  the  series  corresponding  to  this  expression 
will  have  no  terms  in  it,  i.e.,  the  result  is  zero. 

Regardless  of  which  options  are  chosen,  the  final 
solution  of  X,  X,  and  *(=i;)  are  output  to  disk  files  named 
F0R21.DAT,  FOR22.DAT,  and  FOR23.DAT.  This  output  is  in 
coded  form,  i.e.,  for  each  term  in  the  Poisson  series 
corresponding  to  the  solution,  one  line  containing  six 
words  as  described  earlier  is  written  to  the  file.  (Pre- 
ceding the  six  words  in  each  line  is  a nimber  which  is  of 
significance  only  to  the  algebraic  routine.)  The  variables 
appearing  in  the  series  and  their  names  are  shown  at  the 
beginning  of  each  file.  Aside  from  this  list  of  variables 
and  names,  the  file  carries  no  other  information  which 
identifies  the  series.  It  is  the  responsibility  of  the 
user  to  supply  this  identification  and  to  rename  the  files. 
The  contents  of  these  files  may  be  transformed  into  a more 
readable  form  using  the  utility  program  TYPSER,  which  is 
illustrated  in  an  example  below. 
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D.  Warnings  and  Special  Considerations 


1.  While  LNDSTD  does  not  limit  the  number  of  terms 
in  f(x) , the  maximum  degree  in  x of  f(x)  is  32/i  where  i 
is  the  desired  order.  Practical  considerations  further 
restrict  the  length  and  degree  of  f(x).  For  example,  a 

3 

sine  solution  corresponding  to  f(x)  = x-2x  is  35  terms 

long  to  4th  order  in  and  required  about  one  minute  of 

CPU  time.  A sine  and  cosine  solution  corresponding  to  f(x) 

2 3 

f(x)  = ax  + 3x  + bx  , where  a and  b are  constants,  is 
about  700  terms  long  to  4th  order  and  requires  ten  minutes 
of  CPU  time.  In  the  latter  case,  not  only  is  the  execution 
time  long,  but  the  storage  requirements  approach  the  limits 
of  the  DEC  System  10  computer.  It  is  strongly  recommended 
that  the  solution  for  any  f(x)  be  first  run  only  to  1st 
or  2nd  order  to  get  some  idea  of  the  size  of  the  solution. 
For  the  second  case  above,  a 2nd  order  solution  is  less 
than  100  terms  long. 

2.  The  problem  most  likely  to  be  encountered  by 
the  user  is  that  of  insufficient  core  storage  for  the  job. 

If  this  occurs,  LNDSTD  will  type  out  a message  stating  that 
available  storage  has  been  exceeded,  and  execution  will 
terminate.  At  this  point,  disk  files  F0R21.DAT,  FOR22.DAT, 
and  FOR23.DAT  will  contain  the  solution  to  the  highest 
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order  which  was  complete  when  execution  terminated.  In 
order  to  get  a solution  to  higher  order  than  this,  the 
dimensions  in  LNDSTD  must  be  increased  and  the  job  reloaded. 

3.  If  the  0th  order  solution  of  x is  chosen  to 
be  of  the  form  x^  cos  9 then  the  initial  value  of  x is 
zero  and  does  not  appear  in  the  solution.  In  this  case, 
variable  6 is  available  to  the  user,  thus  increasing  the 
number  of  free  variables  to  6,  i.e.,  variables  6-11.  In 
any  case,  the  name  used  for  variable  6 will  be  that 
entered  for  initial  x at  Step  3. 

4.  In  the  present  version  of  LNDSTD,  all  numeric 
constants  must  be  approximated  by  a rational  fraction  N/D. 
Thus  the  largest  constant  possible  is  + 34359738367  and 
the  smallest  possible  is  + 1/34359738367.  If  during  the 
procedure  a number  exceeds  these  bounds,  execution  will 
terminate  and  the  message  "LOSING  ACCURACY  IN  REFACT"  will 
be  typed. 

5.  LNDSTD  does  not  check  to  insure  that  the 
exponents  e^^  and  trigonometric  coefficients  do  not  exceed 
their  bounds  (-32  to  +31  or  -16  to  +15).  It  is  the  respon- 
sibility of  the  user  to  insure  that  this  does  not  occur. 

If  it  does  occur,  the  result  is  that  the  exponent  or  co- 
efficient overflows  into  the  position  of  an  adjacent 
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variable,  and  the  end  result  will  be  incorrect.  It  will  be 
evident  in  the  final  check  if  this  has  occurred.  If  the 
check  is  suppressed  by  choosing  option  1,  then  the  fact 
than  an  error  has  occurred  will  be  evident  only  by  the 
appearance  of  undesired  variables  or  in  the  drastic  effects 
on  the  variable  which  overflowed. 

E . Example  Runs 

Use  of  Programs  LNDSTD  and  TYPSER  are  illustrated 

in  the  three  examples  below.  The  first  example  solves  the 

2 3 

equation  u + u + €(au  + 3u  + bu  ) = 0 to  1st  order  in  €. 

In  the  second  example,  the  disk  file  EXAMPL  is  created  for 

input  to  LNDSTD  to  solve  the  differential  equation 
3 

“+T+e(-n-2n)=0.  In  this  case,  no  options  are 
chosen  and  output  is  to  disk  file  EXAMPL. LP,  which  is  later 
queued  to  the  line  printer.  (Note  that  in  both  examples, 
unused  variables  are  given  the  name,  "????".  Since  the 
appearance  of  these  variables  in  any  of  the  results  indi- 
cates an  error,  this  will  serve  to  flag  such  an  occurrence.) 

Example  3 illustrates  the  use  of  TYPSER.  The  program 
is  self-prompting  and  needs  little  explanation.  If  vari- 
able names  appear  at  the  beginning  of  the  input  file, 

TYPSER  will  use  them.  If  the  user  wishes  to  use  different 
names,  all  lines  containing  alphabetic  characters  must  be 
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removed  from  the  input  file.  TYPSER  will  then  ask  for  new 
names  and  corresponding  variable  numbers,  first  for  the 
polynomial  variables,  and  then  for  the  trigonometric  vari- 
ables. A carriage  return  indicates  that  no  more  names  are 
to  be  entered.  TYPSER  will  also  truncate  a series  to  any 
order  in  any  polynomial  variable  and  will  reverse  the  print 
order  if  desired.  (Series  are  normally  printed  in  order 
of  decreasing  powers  of  variable  1).  The  series  printed 
in  example  3 are  taken  from  examples  1 and  2.  The  variable 
names  in  the  last  series  are  changed  as  illustrated. 
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Example  1: 


.RUN  DSKD:LNDSTD[ 363, 1056] 
LNDSTD:VERSI0N6 REVISED  lIKilll 


START  OF  EXECUTION 


l.IS  DATA  TO  BE  ENTERED  FROM  A TERMINAL(T)  OR  FROM  DISK(D)?  T 
IS  PRINTED  OUTPUT  TO  BE  TO  TERMINAL (T)  OR  TO  DISK(D)?  T 


2.  ENTER  MAXIMUM  ORDER  IN  E OF  FINAL  SOLUTION:  1 

3.  ENTER  NUMBER  OF  TERMS  IN  F(X):  3 

4.  ENTER  DESIRED  NAMES  FOR  VARIABLES  AND  PARAMETERS: 

E-  EPS 
X-  U 
INITIAL  X-  UI 
INITIAL  XDOT-  UDI 
TAU«  TAU 
PARAMl-???? 

PARAM2-  A 
PARAM3-  B 
PARAM4*???? 

PARAM5-???? 


5. ENTER  F(X)  ONE  TERM 
INUTl,  IDNOM,  lEXPl, 
NNNNNNNNNNDDDDDDDDDD 
1 1 

3 1 

1 1 


AT  A TIME.  THE  ELEMENTS  OF  EACH  TERM  ARE 
IEXP2  AND  THE  FORMAT  IS  (2110,2(3X012)): 


000000000000 

204040404040 

204040404040 

204040404040 


000000000000 

204140404041 

204040404042 

204041404043 


6. SINE,  COSINE,  OR  SINE  & COSINE  S0LUTI0N(S,C,  OR  SC)?  SC 

7. ENTER  OPTIONS  (TYPE  "2"  FOR  LIST);  F0P.MAT(61):  2 
OPTl-1  MEANS  SKIP  CHECK  AT  END 
0PT2-1  MEANS  DO  NOT  PRINT  FINAL  X 
0PT3-1  MEANS  DO  NOT  PRINT  FINAL  XDOT 
0PT4-1  MEANS  DO  NOT  PRINT  FINAL  TAUDOT 
0PT5-1  MEANS  DO  NOT  PRINT  INTERMEDIATE  RESULTS 
0PT6-1  MEANS  DECIDE  PRINTING  ON  INDIVIDUAL  BASIS 
ENTER  OPTIONS:  000001 
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HE  PERTURBING  FUNCTION  IS: 


F(X)  - 
+A*U 
+B*U**3 
+3*U**2 


THE  SOLUTION  WILL  BE  A SINE  & COSINE  SOLUTION  TO  ORDER  1 IN  EPS. 
THE  OPTIONS  CHOSEN  ARE: 

OPT6:  DECIDE  PRINTING  ON  INDIVIDUAL  BASIS 

PROCEED (Y/N)?  Y 

SOLUTIONS  FOR  U,  UDOT  AND  TAUDOT  WILL  BE  OUTPUT  IN  PACKED 
FORMAT  TO  FILES  NAMED  FOR21.DAT,  FOR22.DAT,  AND  FOR23.DAT. 


*****************ic*iii,***i,i,*-)t*ii**itie****irk********i€****ii************iiitit*** 

*************************  ZEROTH  ORDER  *************************** 
************************************************************************ 


U0  - ( 2 TERMS) 

PRINT (Y/N) ?Y 
+UI*COS(TAU) 
+UDI*SIN(TAU) 


************************************************************************ 
**************************  FIRST  ORDER  *************************** 

**********************************A************************************* 

OMEGAl  - ( 3 TERMS) 

PRINT (Y/N)?Y 
+3/8*UDI**2*B 
+1/2*A 
+3/8*B*UI**2 


U1  - ( 17  TERMS) 

PRINT(Y/N)?N 

********************************* A4t************************************* 

**************************  FINAL  RESULT  ************************** 
************************************************************************ 
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U =•  ( 19  TERMS) 

PRINT(Y/N)?Y 

-3/32*EPS*UDI**2*B*UI*COS(3TAU) 

■M/32*EPS*B*U1**3*C0S(3TAU) 

-1/32*EPS*UDI**3*B*SIN(3TAU) 

+3/32*EPS*UDI*B*UI**2*SIN(3TAU) 

-l/2*EPS*UDI**2*COS(2TAU) 

+l/2*EPS*UI**2*COS (2TAU) 

+EPS*UDI*UI*SIN(2TAU) 

+3/ 32*EPS*UDI**2*B*UI*COS (TAU) 
+2*EPS*UDI**2*COS(TAU) 

-1/ 32*EPS*B*UI**3*COS (TAU) 
+EPS*UI**2*C0S (TAU) 

+UI*COS(TAU) 

-9/32*EPS*UDI**3*B*SIN(TAU) 

-1/2*EPS*UDI*A*SIN(TAU) 

-21/32*EPS*UDI*B*UI**2*SIN(TAU) 

-2*EPS*UDI*UI*SIN(TAU) 

+UDI*SIN(TAU) 

-3/2*EPS*UDI**2 

-3/2*EPS*UI**2 


UDOT  - ( 17  TERMS 

PRINT (Y/N)?N 


TAU  DOT  - ( 4 TERMS) 

PRINT (Y/N)?Y 
+3/8*EPS*UDI**2*B 
+1/2*EPS*A 
+3/8*EPS*B*UI**2 
+1 


CHECK  *************************** 

************************************************************************ 

D2(X)/DT2  + X + E*F(X;E)  - ( 0 TERMS) 

STOP 

END  OF  EXECUTION 

CPU  TIME:  17.06  ELAPSED  TIME:  5:16.10 
EXIT 

.RENAME  U-FOR21.DAT,  UDOT-FOR22.DAT,  TAUDOT-FOR23.DAT 
Files  renamed: 

FOR21.DAT 

FOR22.DAT 

FOR23.DAT 
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.SOS  EXAMPL 
Input:  EXAMPL 
00100  4 

00200  2 

00300  EPS 

00400  4H 

00500  ???? 

00600  HDI 

00700  TAU 
00800  ???? 

00900  ???? 

01000  ???? 

01100  ???? 

01200  ???? 

01300  1 1 204040404040  204040404041 

01400  -2  1 204040404040  204040404043 

01500  S 
01600  0$ 

*ES 

[DSKD: EXAMPL] 

.RUN  DSKD :L.NDSTD[ 363. 1056] 

LiroSTD: VERSION  6 REVISED  2/16/77 


START  OF  EXECUTION 


l.IS  DATA  TO  BE  EITTERED  FROM  A TERMINAL (T)  OR  FROM  DISK(D)?  D 
IS  PRINTED  OUTPUT  TO  BE  TO  TERMINAL (T)  OR  TO  DISK(D)?  D 
ENTER  INPUT  FILE  NAME:  EXAMPL. 

ENTER  OUTPUT  FILE  NAME:  EXAMPL. LP 
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THE  PERTURBING  FUNCTION  IS: 


F(X)  = 
-2*H**3 
+H 


THE  SOLUTION  WILL  BE  A SINE  SOLUTION  TO  ORDER  4 IN  EPS. 

NO  OPTIONS  WERE  CHOSEN. 

PROCEED (Y/N)?  Y 

SOLUTIONS  FOR  H,  HDOT  AND  TAUDOT  WILL  BE  OUTPUT  IN  PACKED 
FORMAT  TO  FILES  NAMED  FOR21.DAT,  FOR22.DAT,  AND  FOR23.DAT. 


************************************************************************ 

*************************  ZEROTH  ORDER  *************************** 
************************************************************************ 


H0  = ( 1 TERMS) 


************************************************************************ 

*************************  first  ORDER  *************************** 
************************************************************************ 


OMEGAl  - ( 2 TERMS) 

HI  » 3 TERMS) 


************************************************************************ 
*************************  SECOND  ORDER  *************************** 

***************************************H[***********A*****tl,A:«,  ******  ****** 

OMEGA2  - ( 3 TERMS) 


H2  ( 


6 TERMS) 


********:*?******  A **ye******y<j':***yc?'f**?JfA******'^c:^  ********  *^***3ir*yf  ******  ***><?** 

*************************  THIRD  ORDER  *************************** 

************************************************************************ 

OMEGA3=  ( 4 TERMS) 

H3  = ( 10  TERMS) 


************************************************************************ 

************************  FOURTH  ORDER  **************************** 

************************************************************************ 

0MEGA4  =•  ( 5 TERMS) 

H4  - ( 15  TERMS) 


************************************************************************ 

**********************  FIXAL  RESULT  ***************************** 

************************************************************************ 

H * ( 35  TERMS) 

HDOT  - ( 31  TERMS) 

TAU  DOT  - ( 15  TERMS) 


************************************************************************ 

*****************************  CHECK  ****************************** 

******************************************************************* **i^** 

D2(X)/DT2  + X + E*F(X;E)  - ( 0 TERMS) 

**  ZERO  ** 

STOP 
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END  OF  EXECL’TION 

CPU  TIME:  1:9.22  ELAPSED  TIME:  4:59.00 

E.XIT 

.RENAME  H-F0R21.DAT,  HDOT-FOR22.DAT,  TAUHDT-FOR23.DAT 
Files  renamed: 

FOR21.DAT 

FOR22.DAT 

FOR23.DAT 

.PRINT  EXAMPL.LP 

Total  of  11  blocks  in  1 file  in  LPT  request 


Example  3 


.SOS  V=U 
Edit:  U 
* FENDS 
00800  END 
*P0:800 


00100 

THE  VARIABLES  USED  ARE 

1 

00200 

POLYNOMIAL 

VARIABLE 

1 

EPS 

00300 

POLYNOMIAL 

VARIABLE 

6 

- 

UDI 

00400 

POLYNOMIAL 

VARIABLE 

8 

M 

A 

00500 

POL’fNOMIAL 

VARIABLE 

9 

» 

B 

00600 

POLYNOMIAL 

VARIABLE 

12 

S 

UI 

00700 

TRIG. 

VARIABLE 

6 

a 

TAU 

00800  END 
*D0:800 

8 Lines  (00100/1:00800)  deleted 
*ES 


24 


(DSKD: V] 

.RUN  DSKD:TYPSER[363,lid56] 
ST.-yiT  OF  EXECUTION 


ENTER  INPUT  FILE  NAME:  LTX)T. 

ENTER  OUTPUT  FILE  NAME:  UDOT.TTY 
ENTER  SERIES  NAME:  U TO  1ST  ORDER 
TRUNCATION? (Y/N) : N 
REVERSE  PRINT  ORDER? (Y/N):  Y 
CONTINUE  OR  QUIT?(C/0)  C 
ENTER  INPUT  FILE  NAME:  H. 

ENTER  OUTPLT  FILE  NAME:  H2.TTY 
ENTER  SERIES  NAME:  H TO  2ND  ORDER 
TRL’NCATION?(Y/N):  Y 
variable  NL’MBER  and  ORDER:  1 2 
REVERSE  PRINT  ORDER? (Y/N):  N 
CONTINUE  OR  QUIT?(C/Q)  C 
ENTER  INPUT  FILE  NAME:  V. 

ENTER  OUTPLT  FILE  NAME:  V.TTY 
ENTER  SERIES  NAME:  V TO  1ST  ORDER 
INPUT  VARIABLE  NAMES: 

POLYNOMIAL  VARIABLE:  1 EPSV 
POLYNOMIAL  VARUBLE:  6 VDI 
POLYNOMIAL  VARIABLE:  8 ALFA 
POLYNOMIAL  VARIABLE:  9 BETA 
POLYNOMIAL  VARIABLE:  12  VI 
POLYNOMIAL  VARIABLE: 

TRIG.  VARIABLE:  6 TAUV 
TRIG.  VARIABLE: 

TRUNCATION? (Y/N):  N 
REVERSE  PRINT  ORDER? (Y/N) : N 
CONTINUE  OR  QLTT?(C/Q)  Q 
STOP 

END  OF  EXECUTION 

CPU  TIME:  13.73  ELAPSED  TIME:  11:24.23 
EXIT 

.TYPE  UDOT.TTY,  H2.TTY,  V.TTY 
SERIES  NAME:  U TO  1ST  ORDER 
SERIES  LENGTH:  17 
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-UI*SIN(F) 

-EPS*L'I**2*SIN(F) 

-11/32*EPS*B*UI**3*SIN(F) 

-1/2*EPS*A*UI*SIN(F) 

-2*EPS*UDI**2*SIN(F) 

-15/32*EPS*UDI**2*B*UI*SIN(F) 

-t-UDI*COS(F) 

-2*EPS*UDI*UI*C0S(F) 

-9/ 32*EPS*UDI*B*UI**2*COS (F) 

+3/ 32*EPS*UDI**3*B*COS (F) 

-EPS*UI**2*SIN(2F) 

+EPS*UDI**2*SIN(2F) 

+2*EPS*UDI*UI*C0S (2F) 

-3/32*EPS*B*UI**3*SIN(3F) 

■t-9/32*EPS*UDI**2*B*UI*SIN(3F) 

+9/32*EPS*UDI*B*UI**2*COS(3F) 

-3/32*EPS*UDI**3*B*COS(3F) 

SERIES  NAME:  4 TO  2ND  ORDER 

SERIES  LENGTH:  10 

+1/256*EPS**2*HDI**5*SIN(5F) 

+3/16*EPS**2*HDI**5*SIN(3F) 

-5/32*EPS**2*HDI**3*SIN(3F) 

+1/16*EPS*HDI**3*SIN(3F) 

+271/256*EPS**2*HDI**5*SIN(F) 

-45/32*EPS**2*HDI**3*SIN(F) 

+3/8*EPS**2*HDI*SIN(F) 

+9/16*EPS*HDI**3*SIN(F) 

-1/2*EPS*HDI*SIN(F) 

+HDI*SIN(F) 

SERIES  NAME:  V TO  1ST  ORDER 

SERIES  LENGTH:  19 

-3/32*EPSV*VDI**2*BETA*VI*COS(3TAUV) 

+l/32*EPSV*BETA*VI**3*COS(3TAL'V) 

-1  / 32*EPSV*VDI**  3*BETA*SIN  ( 3TAL”V) 

+3/32*EPSV*VDI*BETA*VI**2*SIN(3TAUV) 

-l/2*EPSV*VDI**2*COS(2TAUV) 

+1 / 2*  EPS V*VI**  2*CO  S ( 2TAUV ) 
+EPSV*VDI*VI*SIN(2TAUV) 

+3/ 32*EPSV*VDI**2*BETA*VI*COS (TAUV) 
+2*EPSV*VDI**2*COS (TAUV) 

-1/ 32*EPSV*BETA*VI**  3*COS (TAUV) 
+EPSV*VI**2*COS (TAUV) 

+VI*COS(TAUV) 

-9/32’»EPSV*VDI**3*BETA*SIN(TAUV) 

-1/2*EPSV*VDI*ALFA*SIN(TAUV) 

-21/32*EPSV*VDI*BETA*VI**2*SIN(TAUV) 

-2*EPSV*VDI*VI*SIN(TAUV) 

+VDI*SIN(TAUV) 

-3/2*EPSV*VDI**2 

-3/2*EPSV*VI**2 
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